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ABSTRACT 



Aims. We study the disc planet interactions of low-mass protoplanets embedded in a circumstellar disc. We extend the standard theory 
of planet migration from the usual locally isothermal assumption to include non-barotropic effects, focusing on the validity of linear 
theory. 

Methods. We compared solutions of the linear equations with results from non-linear hydrodynamic simulations, where in both cases 
we adopted a background entropy gradient and solved the energy equation. 

Results. We show that the migration behavior of embedded planets depends critically on the background radial entropy gradient in 
the disc. The presence of such a gradient not only changes the corotation torque on the planet, but also always guarantees a departure 
from linear behavior, which gives a singular density response at corotation, in the absence of thermal or viscous diffusion. A negative 
entropy gradient tends to give rise to positive, non-linear corotation torques apparently produced as material executes horseshoe turns 
at approximately constant entropy. These torques have no counterpart in linear theory, but can be strong enough to push the planet 
outwards until saturation starts to occur after a horseshoe libration period. Increased thermal diffusion acts to reduce these non-linear 
torques, but, at the same time, it can help to prevent their saturation. In combination with a small kinematic viscosity that is able to 
maintain a smooth density profile the positive torque could be sustained. 

Key words. Hydrodynamics - Planets and satellites: formation 



1. Introduction 

, During their formation process inside circumstellar discs, plan- 
I ets can change their orbital pa rameters by gravitational in - 
teraction with the gaseous disc jGoldreich & Tremainelll980l) . 
Torques generated at various resonances can promote or damp 
eccentricity growth jGoldreich & Sarj 2003h. an d change the 
semi-major axis of the forming planet dWardl 19*9 7^. One can dis- 
tinguish torques due to waves, generated at Lindblad resonances, 
which propagate away from the planet, and corotation torques, 
generated near the orbit of th e planet, where material on av erage 
corotates with the planet (see lGoldreich & Tremaindl 19791) . The 
usual result of these torques is that the planet migrates inwards, 
towards the central star, at a rate t hat is determ ined by the planet 
mass and the disc parameters (see lWardlll997h . 

In terms of planet and disc mass, three types of migration 
can be distinguished. Low-mass planets, whose gravitational in- 
fluence is not strong enough to overcome pressure effects, gen- 
erate a linear disc response that gives rise to an inward migration 
rate that is pro portiona l to the planets mass. This is called Type 
I migration (Ward 1997). For standard Solar nebula parameters. 
Type I migration applies to planets up to several Earth masses 
(Me). 

Higher-mass planets excite non-linear waves, and are able 
to tid ally truncate the disc, opening an annular gap around their 
orbit dLin & Papaloizoull979l[T986h . Migration occurring when 
such a gap is maintained is referred to as Type II migration dWardl 
[l997), and it takes the planet inward on a viscous time scale. 
The minimum mass for opening a gap, again for standard disc 
parameters, is approximately one Jupiter mass. 



Both types of migration have been studied extensively 
through numerical hydrodynamical simulations, in a tw o- 
dimensional set-up dNelson et alj2000l:lD'A ngelo et al.'2002') as 
well as in three dimensions ( D'Angelo et al. 2003b; Bate et al] 
'2003), giving good agreement in the respective mass regimes. In 
the intermediate mass regime, depending on the disc viscosity 
and surface density distribution, the corotation torque may get a 
non-linear b oost to the extent t hat it determines the sign of the 
total torque dMasset et alj|2006l) . 

Corotation torques are also responsible for a third 
type of migration, sometimes called runaway migration 
( Masset & Papaloizou 2003). This Type III migration 
dArtvmowicj |2004|) can be very fast, and applies to plan- 
ets t hat open up a parti al gap and are embedded in a massive 
disc dMasset & Papaloiz ou 2003). 

In all studies referred to above, a simplified equation of state, 
for which the pressure depends on the density and the local ra- 
dius only was used. The dependence on radius comes about from 
adopting a fixed radial temperature profile. As at any location, 
the pressure depends only on the density, we describe such an 
equation of state as being locally barotropic. The use of such an 
equation of state removes the need to solve the energy equation 
and thus makes the equations more tractable to tackle both ana- 
lytically and numerically. The question whether this approach is 
valid has only been addressed fairly recently th rough numerical 
simulations that do iiiclude the energy equation dP'Angelo et al.l 
'2003a;'Klahr & Kley 2006"), focusing on high -mass planets only. 
^Morohoshi & Tanaka (2003) studied the effect of optically thin 
cooling on disc-planet interactio ns using a local approach, while 
iJang-Condell & Sasselovl d2005h calculated the torque on low- 
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mass planets analytically in discs with a realistic temperature 
pro file. 

iPaardekooper & Mellemal ( l2006a f) showed that relaxing the 
barotropic assumption can change the migration behavior of 
low-mass planets dramatically. There it was shown, using three- 
dimensional radiation-hydrodynamical simulations, that, for the 
disc parameters adopted, low-mass planets migrate outwards 
through the action of a strong positive corotation torque. Due 
to the excessive computational overheads associated with such 
simulations, it was not feasible to do a parameter survey for dif- 
ferent planet masses and disc models. 

In this paper, we aim to clarify the origin of the positive coro- 
tation torque, its dependence on planet mass and disc parameters 
together with the relationship of this type of disc planet interac- 
tion to the linear perturbation theory that has been successfully 
applied to the Type 1 migration regime. 

To ease the computational demands we focus on non- 
barotropic two-dimensional discs (commonly termed cylindrical 
discs, see iHawlevI (l2000h ) which incorporate thermal diffusion 
rather than the more complex radiative transfer 

The plan of the paper is as follows. In Sect. |2] we review the 
governing equations and the disc models used. Section [3] is de- 
voted to the linear theory of the interaction of a low-mass planet 
with a gaseous disc, and we solve the resulting equations in Sect. 
|4] In order to remove singularities we adopted the well known 
Landau prescription together with a small thermal diffusivity ap- 
plied in some cases. Although migration torques could be calcu- 
lated in the limit of zero dissipative effects, the density response 
became singular implying that non-linear effects may always be 
important in the corotation region. 

After describing the numerical set-up in Sect. |5j we perform 
fully non-linear hydrodynamical simulations in Sect. |6] Indeed 
we find that a non linear small scale coorbital density structure is 
set up at early times. We go on to investigate the effect of thermal 
diffusivity, finding that a non zero value is required to ensure the 
coorbital region is adequately resolved. The associated torque is 
then indeed found to be positive, for favorable density gradients, 
and sufficiently strong and negative entropy gradients, for about 
one horseshoe libration period after which saturation sets in. In 
Sect. |2] we present a simple non-linear model of the corotation 
region based on material executing horseshoe turns with con- 
stant entropy, which can result in a density structure that leads to 
a temporary positive torque on the planet when there is a back- 
ground negative entropy gradient and compare it with simula- 
tions. 

Whether the positive torque can be sustained for longer times 
will depend on the evolution of the planetary orbit and disc and 
whether these circumstances provide an appropriate corotational 
flow that can resupply appropriate low entropy material. To in- 
vestigate this aspect, we study the long-term evolution of the 
torque in Sect. [8] and in the case of a particular example, show 
that for an appropriate rate of thermal diffusion and a small vis- 
cously driven mass flow through corotation, the entropy-related 
corotation torque can be sustained for several libration periods. 

We go on to give a discussion of our findings in Sect. |9]and 
conclude in Sect. [TO] 

2. Basic equations and disc modeis 

The basic equations are those of the conservation of mass, mo- 
mentum and energy for a two dimensional cylindrical disc in a 
frame rotating with angular velocity ftp in the form 

dp 
di 



— + 2f}pk X V 



DE 



-VP- V$ 
P 



P Dp DS 

p - = pT = -V • F. 

Dt p Dt ^ Dt 



(2) 
(3) 



Here, p denotes the density, v = {vr,v,p) the velocity, k de- 
notes the unit vector in the vertical direction, P the pressure and 
$ — $G + f^p?'^/2 where $g is the gravitational potential. The 
convective derivative is defined by 

D d 

and F = —KWT is the heat flux, with K being the thermal con- 
ductivity and T being the temperature (of course radiation trans- 
port may be incorporated in this formalism). In our numerical 
simulations, we choose K = K{r) such that the initial temper- 
ature profile gives V • F = 0. We adopt an ideal gas equation of 
state such that 



P = RgpT, 



(5) 



with Rg being the gas constant. The internal energy per unit mass 
is given by 



E 



P 



(7-l)p' 



(6) 



where 7 is the constant adiabatic exponent and 

= Rg \og{P/p'') /(7 — 1) is the entropy per unit mass. 
We adopt a cylindrical polar coordinate system (r, (p) with 
origin (r = 0) located at the central mass. The self-gravity of 
the disc is neglected. Thus the gravitational potential is assumed 
to be due to the central mass and perturbing planet such that 
$G *G0 + *Gp, where 



$G0 = 



and 



$Gp = 



GIVL 



(7) 



-GAL 



r'^ + r'i — 2r7'p cos{(p — (pp) + b'^r. 



+ 



GMpT C0s((y9 — (pp) 



(8) 



-V • (pv) 



(1) 



In the above the last term is the indirect term, denotes the 
mass of the central object, with Mp, Tp, and ipp denoting the 
mass, radial and angular coordinate of the protoplanet respec- 
tively. We also allow for a gravitational softening parameter h. 

We remark that there are several limiting cases in the above 
description. When heat transport is neglected, {K = 0), the sys- 
tem is adiabatic. When the energy equation is dropped and the 
temperature is taken to be a fixed function of r we have the usual 
locally isothermal limit. When fip = the reference frame is 
non-rotating but non-inertial as the origin accelerates together 
with the central mass. This is accounted for by the indirect term 
in the potential. Numerical calculations are most conveniently 
performed in a frame corotating with the protoplanet. Then Vlp 
becomes the circular Keplerian angular velocity at radius Tp. At 
a general radius, r, this is given by fix = (GAf^/r'^)^/^. 

The formulation given above applies to a cylindrical disc 
model where there is no vertical stratification or dependence on 
the vertical coordinate. The equations apply to a cylindrical sys- 
tem with a constant vertical semi-thickness, iJo which may be 
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Fig. 1. The functions U and Q divided by the protoplanet to central star mass ratio are plotted for po oc r~^/^, b = 0.03 and m = 2 
in dimensionless units. Real parts are plotted with dashed lines , imaginary parts with full lines. Curves for the Landau parameters 
e = 10~^ e = 3. X 10~^ are shown. The response is very similar for these cases but the latter case is slightly more damped which 
is manifest by somewhat smaller amplitude swings. 



specified arbitrarily as there is no explicit dependence on it. The 
density, p, is taken to have the same spatial variation as the sur- 
face density to which it is related through E = 2pHo. For a 
thin disc of the type we wish to consider, at any location there 
is a putative local vertical semi-thickness H — Cs/^k, with 
Cs = \fPjp being the local isothermal sound speed that is a 
function of position and associated with the neglected vertical 
stratification. When considering the physical state variables as- 
sociated with any such location we regard as having been 
adjusted to coincide with H. 



2. 1. Equilibrium models 

In this paper we work with equilibrium model discs into which a 
perturbing protoplanet is subsequently inserted. These are such 
that density has a power law dependence of the form p cx r^", 
with a being constant. We also adopt T cx r^^, giving a pu- 
tative semi-thickness H (x r, then we have P oc r^"^^. The 
local gas angular velocity required for radial hydrostatic equilib- 
rium is given hy il — ilK(l — {a + l){H/r)'^). Thus for these 
models fl cx f^K and the epicyclic frequency k — fl. Here we 
recall in general that = (2r2/r)(i(r^ri)/(ir. We also adopt the 
representative value 7 = 1.4 unless stated otherwise. 



3. Linear theory 

We here consider the response of the disc to the forcing of a 
protoplanet in a circular orbit assuming that the induced pertur- 
bations are small so that the basic equations may be linearized. 



We perform a Fourier decomposition of the perturbing potential 
such that in the non-rotating frame 




Gpm exp {inup — imflKt) 



(9) 



where TZe denotes that the real part is to be taken. Here, as well 
as in the rest of this section, JIk is short for ^K{rp)- 

Similarly, the perturbation response of state variables that 
are non-zero in the equilibrium, indicated with a prime, for each 
value of TO, has an exponential factor exp {irmp — imflKt) with 
an amplitude depending only on r. Equations for these ampli- 
tudes are obtained by linearizing the basic equations. We begin 
by considering the adiabatic limit in which the thermal diffusiv- 
ity is set to zero. Then we linearize the adiabatic condition 



DS dS „ 
to obtain 

5' - ^r- dS 

ia dr ' 



(10) 



(11) 



where cr = m{il — JIk)- Expressed in terms of the density and 
pressure perturbations Eq. (fTTT l takes the equivalent form 



p_ 

p 



VrA 

iagr 



(12) 



where gr = — ( 1 / p) {dP/ dr) and the square of the Brunt- Vaisala 
frequency is given by 



^ _ IdP / 1 dP I dp 
p dr \jP dr p dr 



(13) 
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Fig. 2. The density response divided by the protoplanet to central star mass ratio, for po oc r^/^ , 6 = 0.03 and to = 2 in dimension- 
less units. Real part dashed line , imaginary part full line. For the left panel the Landau parameter e = 10^^ and for the right panel 
e = 6 X 10~^. In the former case the amplitudes of the real and imaginary parts reach extreme values of order 10^. 



which can be negative in regions of formal local convective in- 
stability. Equation (fTZt can be used to express the density per- 
turbation in terms of the pressure and velocity perturbations. 

Linearization of the two components of the equation of mo- 
tion together with the continuity equation then yield, after elimi- 
nation of the azimuthal velocity perturbation, a pair of first order 
ordinary differential equations for the quantities Q = P' I P^h 
and U — rP^/'^Vj. which may be written in the form 



— = CiVr 

dr 



dU_ 

dr 



DlQ + D2Ur + S2, 



(14) 



(15) 



where the coefficients are given by 

Ci = -ipP-'^^^^\a-{K^ +A)/a) 
C2 = -~{2mn)/{ra) 
Di = ~ir{a^ - -^Pm^/{r^p))/{a-fP^-^'^) 
D2 = {P^/^mK^)/{2na) 
Si - -pP(-i/^)[(d$Gpm/dr) + (2TOr!)$Gpm/(m)] 
^2 = z(pi/^m2)/(m)$Gpm. 

3. 1. The corotation singularity 

The perturbation of the protoplanet causes the excitation of out- 
going density waves that are associated with a conserved wave 
action or angular momentum flux. This causes a torque to act 
on the protoplanet. However, angular momentum exchange be- 
tween protoplanet and disc may also occur directly at corotation 
where 0. In linear theory, this type of exchange is localized 



at corotation through the operation of a corotation singularity. 
In order to study the domain near corotation where ct = 0, we 
neglect in the first set of brackets in the expression for Di 
above. Then we may derive a single equation for U from Eqs. 
(O and ( fTSb in the form 

P^/t d / pr fdU 
pr dr \ p2/7 dr 



U 



rnP^d_ / pn^ \ 
pra dr \ 2VlP'^h ) 

argr 



(16) 



We see that Eq. JTSI l is in general singular at corotation with 
a second order pole at = 0. This will result in angular momen- 
tum exchange with the perturber In order to be singularity free 
we require both that the entropy gradient be zero or equivalently 
A^Q, and that d{pK.'^ /{2QP^/'^)/dr = 0. The latter condition 
is the generalization of the condition that the gradient of specific 
vorticity, or vortensity, given by 



(17) 



dr \2np 



should be zero that applies in the strictly barotropic case. 



3.2. The Landau prescription 

In order to solve the forced response problem, we have to deal 
with the corotation singularity. Physically this is resolved by in- 
cluding the dissipative effects of heat transport and viscosity. 
However, by increasing the order of the system of equations, this 
significantly increases the complexity of the problem. For many 
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Fig. 3. The cumulative torque in units of g^S(rp)r2^rp acting on the protoplanet plotted as a function of dimensionless radius for 

po oc r^^/^, b — 0.03 and to = 2. For the left panel plots are given for the Landau parameters e — 10~^ (dashed line) e = 1G~^ 
(dotted line) e = 3 x 10^^ (full line). For the right panel plots are given for the dimensionless thermal diffusivities — 10^^ 
(dotted line) and Dc — 10^^ (full line). The cumulative torque is defined to be zero at the outer boundary and equal the total torque 
at the inner boundary. In all of these cases the total torque is positive. 



purposes this can be avoided by using the Landau prescription. 
According to this the forcing frequency is given a small positive 
imaginary part such that a — > m{i^ — JIk — iei^K), where e is 
very small. In practice, as is the case with the problem consid- 
ered here, the torques between disc and protoplanet converge as 
e ^ even though the solution becomes singular. 

The latter feature means that even though torques may be 
convergent the response becomes non-linear near corotation 
making the linear theory invalid. The time l/(f7Ke) may be in- 
terpreted as the time scale over which the perturbing potential 
is turned on. Thus a non-linear breakdown for e below a certain 
value, eo , means that, viewed from the point of view of an initial 
value problem, the linear solution should be considered as being 
vaUd only for a finite time period, which may be estimated as to 
be of order 1/{^Ik^q). 

For the problem on hand, non-linear breakdown occurs 
mainly through terms associated with the entropy gradient or 
A, these being associated with a second order pole. This type 
of singularity may be removed by introducing heat diffusion so 
that we would expect that for a sufficiently high thermal diffu- 
sivity, the validity of linear theory should be restored. This we 
explore below. 



3.3. The effect of a small thermal diffusivity 



In order to investigate the effect of thermal diffusion on the coro- 
tation singularity, we now consider the modification of the con- 
dition given by Eq. ( fTTT ) that occurs when a small thermal diffu- 
sivity is present. 



In this case we must linearize the full energy equation, Eq. 
(|3]l, which leads to 

,^S' + Vr^ = ^V'r. (18) 

ar pi 

Here, because of rapid variation near corotation, when consid- 
ering terms involving K, we neglect everything other than the 
second derivative term. Assuming that the variation of the pres- 
sure perturbation can be neglected, (this approach can be vali- 
dated by inspection of the form of the solutions), we then have 
S' = {dS/dT)pT', and accordingly 

^aS' + Vr^^^V'S', (19) 
ar Cpp 

where Cp = 7Rg/ (7 — 1) is the specific heat at constant pres- 
sure. As we are interested in a local region around corotation, we 
set r — rp + X, where x <C rp, and a — > —SmU.KX / {2rp) . In 
addition we replace {<P /dx"^ — m? /r^). 

Then Eq. (fT9] l can be written in the form 

- i(3zS' = F, (20) 

where z — x + 2Dmi/{3rpilK), P — —3rn^K/{2rpD), 
D = K/{pCp), andF — Vr{dS/dr)/D, the last three quantities 
being assumed constant. Equation ( l20b can be solved in terms 
of in- homogeneous Aiiy functions (e.g.. lAbramowitz & StegunI 
119701) in the form 




Fig. 4. The cumulative torque in units of q'^T,{rp)il.'^rp acting on the protoplanet as a function of m for po oc r 6 = 0.03 (left 
panel) and b = 0.01 (right panel). For both panels plots are given for the Landau parameters e = 10^^ ( full line), e — 10^^ (dotted 
line), e = 10~^ (dashed Une), e = 3 x 10^^ (dot dashed Une). 



with 



c = - 



2{ia + ei){9m/2Dc)^/^ 



(22) 



and 



(23) 

with Dc being the dimensionless diffusivity = -©/(rpfix). 
Here all quantities are evaluated at corotation, r = Vp, apart from 
a SmflKX / {2rp) . 

The in- homogeneous Airy function Hi{C) (see 
lAbramowitz & Stegun 1970) is defined by 



1 f°° 

m(C) = - exp (-fc3 + kQ)dk. 



(24) 



Using this solution, S' is no longer singular at corotation. By 
comparing Eq. ( fTTl i for S' with the solution determined by Eq. 
(I2TI) we can see how to remove the corotation singularity. 

It is apparent that, where it appears multiplied by the entropy 
gradient, the singular quantity 



a 2T:i / 9to 
CT 3m V 2Z?o 



1/3 



Hi 



-2{ia + ei){9m/2Dey/^ 



(25) 



In this way the divergence at ct = associated with terms that 
involve the entropy gradient is removed. When we implemented 
cases with thermal diffusivity as described below, wherever 
otherwise occurred, we applied the Landau prescription with e = 
10"^ 

From Eq. ( l20l i. the length scale associated with the diffu- 
sively controlled corotation region is l/?]^^/"^ ^ rp(£'c/?7i)^^'^. 
Thus for sufficiently large D^, we expect that the effects of coro- 
tation associated with a radial entropy gradient should be re- 
duced such that the validity of linear theory is restored as far 



as these are concerned. Additional restrictions may result from a 
radial vortensity gradient. 



4. Linear response calculations 

We have solved Eqs. ( fT4l i and ( fT5] l for a variety of disc models, 
softening parameters and values of to. We have used both the 
Landau prescription and the combination of that, together with 
a finite thermal diffusivity used to deal with terms associated 
with the entropy gradient, to deal with the corotation singularity. 
The equations were integrated using a fifth order Runge-Kutta 
scheme and outgoing wave conditions were applied. In practice 
these were applied at radii that became closer to the protoplanet 
as TO increased. 

For convenience throughout, we adopt a system of units for 
which Tp = 1, flK{rp) = 1, and the density of the unperturbed 
disc at the protoplanet orbital radius p(rp) = 1. The orbital pe- 
riod of the protoplanet is then 2tt. In all cases the aspect ratio 
Cs/ {rfl) was taken to be 0.05. 

4.1. Results 

The functions U and Q are plotted for the disc with po oc r^^/^, 
b ~ 0.03, being comparable to H/r, and to = 2 in Fig. [T] 
In these cases a Landau prescription was used to deal with the 
corotation singularity. We show curves for the Landau parame- 
ters e ~ 10^^ and e = 3 x 10^^. As is indicated in Fig.[T] the 
response functions are very similar in these cases for which the 
Landau parameter varies by a factor of 3000. The e = 3. x 10^^ 
case not unexpectedly shows slightly more damping. 

However, the density response exhibits different behavior. 
This is plotted for the Landau parameters e = 10~^ and e = 
6 X 10^^ in Fig. |2] for a disc with a = —1/2 (very similar 
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Fig. 5. The cumulative torque in units of g^S(rp)ri^rp acting on the protoplanet as a function of m for po cc r and b = 0.03 

(left panel) and for po oc r~^/^, b — 0.03 but with thermal diffusion (right panel). For the left panel plots are given for the Landau 
parameters e = 10^^ ( full line), e = 10^"* (dotted line), e = 10^'^ (dashed line), e = 3 x lO^'^ (dot dashed line). In the right panel 
plots are given for the dimensionless thermal diffusivities Do = lO^'' (dotted line) and = 10^^ ( full line). The curves almost 
coincide for these cases. 



behavior occurs for a — 1/2). This response shows a strong sin- 
gularity as e is decreased to zero. The extreme values reached 
are c>c 1/e. As the singularity is approached the radial width oc e. 
The effect of this is that although the pressure and velocity re- 
sponse as well as the migration torque converge as e — * 0, the 
density response becomes increasingly singular This means that 
linear theory always breaks down in that limit. This breakdown 
is in fact associated with terms cx A, the entropy gradient and it 
can be removed by considering thermal diffusion. To deal with 
this, one can use the relative smoothness of the pressure and ve- 
locity perturbations to obtain a simple description governed by 
a second order ordinary differential equation (see above). In or- 
der to remain in a regime for which linear theory is valid as far 
as these are concerned it is essential to either incorporate a non 
zero thermal diffusivity or restrict consideration to values of e 
that are high enough for linear theory to be valid. The latter 
viewpoint can be regarded as restricting consideration of evo- 
lution times for an initial value problem to values smaller than 
^ l/(el7K(rp)). 

4.2. Migration torque 

The torque acting on the protoplanet is calculated directly by 
summing the torque integrals evaluated for each Fourier compo- 
nent in the form 



We use the above to define a cumulative torque as a function of 
r for each ni. Thus 



T(rin) = 2HJ„ 



TTmrp'^Gpmdr 



(26) 



where the integral is taken over the unperturbed disc domain 
{Tin, Tout) and Tm denotes that the imaginary part is to be taken. 
The factor 2H takes account of the vertical direction by assum- 
ing the cylindrical disc model to apply over a vertical extent 2H. 



Tp{r)dr, 



(27) 



where rp(7-) is the torque density. Then the total torque acting 
on the protoplanet is T(7-in). 

The cumulative torque acting on the protoplanet is plotted as 
a function of dimensionless radius for po c< r^^/^, 6 = 0.03 and 
TO = 2 in Fig. |3] Plots are given for cases using the Landau pre- 
scription with parameters e = 10^^, e = 10^^ and e = 3x 10^^. 
For comparison plots are given for the dimensionless thermal 
diffusivities Dc = 10^^ (dotted Hne) and Dc — 10^^ imple- 
mented as described above. In all of these cases the net torque is 
positive. We remark that the curves for e = 10~^ and e = 10^* 
are very close as are those for Dc — 10^^ and = 10~^. This 
indicates good convergence of the torques as dissipative effects 
are reduced. 

To emphasize this feature we show the cumulative torque 
acting on the protoplanet as a function of to for po oc r^^/^ 
for the two softening parameters b — 0.03 and b = 0.01 in 
Fig- SI The softening parameter that should be used is uncertain 
but it should be of the order of the aspect ratio to simulate 3D 
effects. Plots are given for the Landau parameters e = 10~^, 
e — 10~^, e — 10^^ and e = 3 x 10^^. The total torques are 
always negative corresponding to inward migration. Again the 
convergence for small e is good but note that torques of larger 
magnitude are obtained for the smaller softening parameter and 
larger m need to be considered. 

Similar results are obtained for different disc models and 
when thermal diffusion is employed as illustrated in Fig. |5] 
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Similarly for a pure Landau prescription, we obtain 




Fig. 6. Contour plot of the full linear density response for po oc 
r^^/^ and thermal diffusivity Dc ~ 10^^ in dimensionless 
units. The linear theory is only likely to be valid for planet 
masses significantly less than an earth mass in this case (see 
Eq. (|28]|). Note the narrow high and low density regions on the 
corotation circle near the protoplanet. These, also seen in the 
non-linear simulations, respectively lead and trail the protolanet 
giving rise to a positive corotation torque. However, the effects 
of this torque are not strong enough to counteract the negative 
Lindblad torques contributed by the high density wakes. 



There we plot cumulative torques as a function of m for po oc 
r~'^/^ and b — 0.03. We show also the corresponding plots for 
Po r"^/'^, b = 0.03 but with the implementation of thermal 
diffusion. 

As indicated above, linear theory always breaks down for 
sufficiently small dissipation parameters, and it is important to 
estimate parameter regimes where it could be valid. We have 
done this both when a pure Landau prescription is used and also 
when thermal diffusion is applied. To do this we calculate the 
full response by summing over m and determine the condition 
that the perturbation to the entropy gradient be less than or equal 
to the same magnitude as the unperturbed value. Clearly non- 
linear effects may set in for weaker perturbations than those that 
result in this condition being marginal. Thus we estimate that 
non-linearity sets in for smaller diffusion coefficients than those 
obtained when this condition is marginally satisfied. We show a 
contour plot of the full linear density response for po cx r^^/^ 
and thermal diffusivity ~ 10^^ in Fig.|6] From this and other 
similar responses we obtain the condition 



'0.03 



1/3 



(L55X (g/10-^)) 
(I?e/10-6)2/3 



0.03 



1/4 



2.5 X Wq 



< 1. 



(29) 



(28) 



Interestingly, from Eq. ( I28l l we estimate that for protoplanets in 
the earth mass range, if Dc < 10~^, there should be departures 
from linear theory. Similarly from Eq. ( |29l ). we estimate pro- 
toplanets in this mass range to be in the non-linear regime if 
e <~ 10^^. This indicates that linear theory may only be rel- 
evant for short evolutionary times in cases with very low dissi- 
pation and then from Fig.|4]we would expect that the full linear 
torque is never established. Although the above estimates are 
uncertain, we comment that the scaling implied by Eq. ( |29] l. that 
for the same degree of nonlinearity e cx 9^/^, implies that the 
time required to attain the same degree of departure from linear 
theory should scale as q^^/^. 



5. Numerical hydrodynamical simulations 

In this Section, we describe the set-up for our numerical hydro- 
dynamical simulations of two-dimensional gas discs with em- 
bedded planets, using a non-barotropic equation of state. The 
initial disc models and system of units are the same as those 
employed for the linear calculations. The softening parameter 
is 6 = 0.03. This approximately corresponds to the putative 
disc semi-thickness and the use of such a softening parameter 
approximates the result of appropriate vertical averaging of the 
gravitational potential. There is no explicit viscosity in the model 
unless otherwise stated. In Sect. l6.3l we include explicit heat dif- 
fusion. For locally isothermal simulations, we drop the energy 
equation and instead use a fixed temperature profile that gives 
rise to a constant aspect ratio in the initial state. For models in- 
cluding the energy equation we have used different temperature 
profiles to vary the entropy gradient while keeping the density 
gradient constant. Unless stated otherwise, we slowly build up 
the mass of the planet during the first three orbits to avoid to 
avoid transients due to the sudden introduction of the planet into 
the disc. For low-mass planets in isothermal simulations this is 
usually unnecessary, but the non-linearities in the flow associ- 
ated with an entropy gradient can introduce some artifacts (see 
below). 



5. 1. Code description 

We use the RODEO method dPaardekooper & Mellem"all2006bh 
to evolve the two-dimensional Euler equations in cylindri- 
cal coordinates (r, (/s). The grid consisted of 512 radial cells, 
equally spaced between r — 0.4 and r — 1.8, and 2048 az- 
imuthal cells, equally spaced over 2tt . We used non-reflective 
boundary conditions (for details see Paardekooper & M ellemal 
2006b) . The energy equation is a straightforward extension of 
the met hod (see [E ulderink & Mellema 1995), which was also 
used in iPaardekooper & Mellemaj (i2006a ) to perform three- 
dimensional, radiation-hydrodynamical simulations of embed- 
ded planets. In this paper, we have not used the radiation dif- 
fusion solver. Heat diffus ion is incorporated the same way as 
Navier-Stokes viscosity in IPaardekooper & Mellemal ( l2006bl) . 

In calculating the torque on the planet, we include all disc 
material, including that within the Hill sphere of the planet. For 
low-mass planets considered here, this should not influence the 
results. 
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Fig. 7. Total torque in units of q^'S{rp)il.^rp as a function of 
time for two different planet masses and two density profiles for 
the disc. The disc was kept locally isothermal, and a viscosity 
1/ = 10~^ ^p^p- The dotted lines indicate the results of our 
linear calculations. 



5.2. Comparison to previous results 

Before embarking on a study of non isentropic discs in Sect. |6l 
we compare our results for locally isothermal discs with previous 
work and expectation from linear theory in this section. 

Different kinds of comparisons of torque calculatio ns with 
linear theory have been reported in the literature. Ma s set et al.l 
(120061) compare locally isother mal three dimensio nal calcula- 
tions with the linear results of iTanaka et al. (20021). This com- 



parison is also reported in lPapaloizou et all 
parison the expression for the Unear torque given by 



2007h . In this com- 



Tl = -(1.35 + 0.55a)q'{H/r)-'j:{rp)ni 



2 „4 



(30) 



is used. On the other hand the results of two dimensional sim- 
ulations using a softened planet gravitational potential, as is 
done here, have also been compared to those given by Eq. 
( [30l l and it has been found that an approximate match may 
be achieved with an app ropriate choice of softening parame- 
ter (Nelson & Papaloizou 2004). We shall show that our results 
are ful ly consistent with those reported in iNelson & Papaloizoul 
dlOOl) and Eq. ^ as well as our own customized linear calcu- 
lations below. We note that the comparisons referred to above 
are done with constant kinematic viscosity coefficients, ly = 
lO^^r-pilp, and so we adopt this value here to make our compar- 
isons. We comment that for this value of v, corotation torques 
are expected to be largely unsaturated. 

We have performed locally isothermal simulations of planets 
with different mass ratio q for different values of a, H/r = 0.05 
and the stated value of ly. 

In Fig.|7]we show the total torque, divided by the mass ratio 
squared, as a function of time for planets with q = 6.34 x 10^^ 
and q = 1.26 x 10^^ for two different values of a. In these cases 
the torques attain almost steady values after about ten orbits with 
any remaining transients due to the initial conditions producing 
only small effects. If the response were fully linear, there would 



be no differences between planets of different masses. It is clear 
that the disc with a = 3/2, which has no vortensity gradient, is 
consistent with the interaction being linear. The torques agree to 
within 10 % with the value obtained from our linear calculations, 
as well as with that obtained from Eq. ( l30l l. For the case with 
a = 1/2, the differences between the two planet masses are 
more pronounced We find the magnitude of the torque divided 
by q^ to be about 10% larger for the smaller mass indicating the 
presence of weak non linear effects. For that mass the torques 
again agree to within 10 % with the value obtained from our 
linear calculations, as well as that found u sing Eq. (l30l l. 

In agreement with Nelson & Papaloizou (2004), torques can 
be fitted by Eq. (|30] l to within 10 — 15% for the mass range con- 
sidered for an appropriate softening parameter. In our case this 
would be close to b = 0.6H/r the value we used for almost all 
cases. [Nelson & Papaloizoul (l2004ft adopted b — 0.7 H/r, which 
is not a significant difference in this context. 



6. The development of positive corotation 

torques resulting from initial entropy gradients 

In this section we focus on the development of positive torques 
arising from initial entropy gradients within the first horseshoe 
libration period, dealing with the possibility of sustaining such 
torques over longer time scales in later sections. We will con- 
sider various disc models, ranging from one with a — —1/2, 
which has strong gradients in entropy and vortensity and a sur- 
face density that increases outwards, to discs with constant en- 
tropy and vortensity and outwardly decreasing surface density. 

6.1. An illustrative case 

We begin by considering a planet with q = 1.26 • 10^^ (cor- 
responding to a 4 planet aroun d a Solar mass star). This 
closely resembles the case studied in Paardekooper & Mellemal 
(l2006al) . The initial density and temperature structure is charac- 
terized hy a = 1/2 and constant H/r. The power law index for 
the entropy is ^ = d{log P/ p'') / dlogr — —0.8, which corre- 
sponds to 7 = 1.4. In Fig.|8]we give contour plots of the density 
(multiplied by r^/^) for an isothermal equation of state and an 
adiabatic simulation after 20 orbits. The adiabatic case can be 
compared with the plot in Fig.|6]obtained from the linear theory 
with some heat diffusion added to regularize the density pertur- 
bation. 

The wave structure differs slightly between the simulated 
adiabatic and isothermal cases. Note that the color scale has been 
chosen to highlight corotation features, but still the differences 
between the spiral waves can be seen. The pitch angle differs 
such that the waves are more open and weaker in the adiabatic 
case. This is due to the fact that linear waves are isentropic, 
leading to a higher sound speed and less compression for equal 
forcing in the adiabatic case. This results in a reduction in the 
strength of the Lindblad or wave torques in adiabatic simulations 
compared to isothermal ones. 

The two simulations in Fig. [8] also differ around corotation 
(near r = rp). The adiabatic simulation shows a strong, nar- 
row density feature which is associated with a positive pertur- 
bation (light shade) for (p > (pp and associated with a negative 
perturbation (dark shade) for ip < <y9p. As the planets are mov- 
ing upwards in Fig. [8j it is clear that this feature acts as to ac- 
celerate the planet in its orbit, causing outward migration. It is 
also found in the linear calculations but at reduced strength (see 
Fig. |6]l. In Sect. Q we will see that this is because linear theory 



10 



S.-J. Paardekooper & J.C.B. Papaloizou: Disc planet interactions in non-barotropic discs 





DSS 


UQD 


1.D2 


1,D4 




Fig. 8. Overview of pr^/^ for a q = 1.26 • 10 ^ planet (4 around a Solar mass star) after 20 orbits. Left panel: locally isothermal 
equation of state. Right panel: adiabatic equation of state. 



fails to take into account the horseshoe turns. We will see below 
that when non-linear simulations are undertaken, this corotation 
torque can easily come to dominate the Lindblad torques, and 
that the sign depends on the direction of the entropy gradient in 
the disc. The isothermal simulation also shows a corotation fea- 
ture, which is due to the radial vortensity gradient in the disc. 
Although the sign of this corotation torque is positive, it is not 
strong enough to overcome the negative wave torque. For plan- 
etsof higher mass the situation may be different (Mas set et al.i 
l2006h . 

In Fig. |9]we compare the time evolution of the total torque 
for an isothermal and an adiabatic simulation with a = —1/2, 
together with a run with the same density profile, but with a 
temperature profile such that the entropy is constant. In the case 
with a negative entropy gradient, we find a positive torque on the 
planet after approximately 10 orbits. After that time, the torque 
continues to rise steadily in all cases. This is when the corotation 
torque is set up, see Fig. [TO] which shows the torque associated 
with the m — 2 component of the density perturbation for the 
case with negative entropy gradient. The torque difference be- 
tween 20 and 10 orbits arises entirely at corotation. In all cases 
shown in Figs. |9] and [TOl the corotation torque is expected to be 
positive due to the positive vortensity gradient and negative (or 
zero) entropy gradient. 

However, in linear theory the magnitude of the corotation 
torque has been found to be not enough to produce a positive 
torque and thus smaller than the wave torque. We note that one 
finds from Tanaka et al. (2002) that the positive linear corotation 
torque for a — —1/2 in the isothermal case is expected to be 
54% of the wave torque in magnitude making the total torque 



46% of the latter. Also, the linear torques are expected to reach 
their final value on a dynamical time scale, which is independent 
of the planet mass ratio, rather than over 20 orbits as is seen 
in Fig. |9] The precise time to reach the maximum depends on 
the time scale on which the planet is introduced, but it is always 
significantly more than a dynamical time scale. Below, we will 
see that this time depends on the planet mass ratio and so this 
behavior must be associated with non-linearities in the flow and 
we interpret it as arising from torques generated after horseshoe 
turns as described in Sect. |2]below. 

After approximately 35 orbits, the torque starts to drop again 
in all cases. This is also an indication that corotation torques play 
a major role, because it is expected that they start to saturate after 
half a libration period, which is close to 35 orbits for this planet. 
We note that the torque then starts to oscillate on the libration 
time scale and, ultimately after the action of non-linear effects 
and diffusion, the disc settles down to a state with zero entropy 
and vortensity gradients around corotation in an average sense. 
Therefore, the corotation torque vanishes unless some form of 
diffusive or evolutionary process can restore the entropy gradient 
within one libration period. We will come back to this issue in 
Sect. [8] 

6.2. Non-linear evolution 

Even at early times the evolution of the torque is dominated 
by non-linear effects. In Fig.[TT]we show the time evolution of 
the torque for three different planet masses. Torques have been 
scaled by q^^, so that if the response was purely linear all curves 
should lie on top of each other The time scale for the rise to- 
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Fig. 9. Total torque in units of q^Y,{rp)fl^r'^ as a function of 
time for a q = 1.26 • 10^^ planet (4 around a Solar mass 
star), using a locally isothermal equation of state (solid line) and 
solving the full energy equation, for two different radial entropy 
gradients ^ ~ d\og{P / p"*) / d\ogr (dotted line: ^ = 0, dashed 
line: ^ = -1.2). All models have a — -1/2. 
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Fig. 10. The cumulative torque in units of q^Yi{r-p)VL\r^ acting 
on the protoplanet plotted as a function of dimensionless radius 
for po oc r^/^, £^ = —1.2 and to = 2 at different times. The 
cumulative torque is defined to be zero at the outer boundary 
and equal the total torque at the inner boundary. 



wards positive values clearly depends on the mass ratio which 
would not be the case if the effect was linear. The time to reach 
the same degree of non-linearity is consistent with being pro- 
portional with q^^/^, as implied by Eq. (|29] l. Similarly, the time 
scale on which saturation sets in, resulting in the torque drop- 
ping below zero again, depends on q. Therefore, although there 
is a linear effect associated with the entropy gradient, the phe- 
nomenon that results in the torque switching sign is clearly non- 
linear. 

This is further illustrated in Fig. [121 where we show the 
TO = 2 component of the surface density perturbation due to 
a g = 3.17 • 10^^ planet (corresponding to a 1 planet 
around a Solar mass star) after 10 and 20 orbits. The imaginary 
parts (solid lines) again are directly related to the torque. At the 
boundaries of the domain we see clear signatures of outgoing 
waves. The difference between t — 2Q and t ~ \Q arises near 
corotation only (see also Fig. [TOli. At t = 10, the perturbation 
is reasonably smooth and resembles the fully developed linear 
case away from corotation (see Fig.|2|. As time progresses, how- 
ever, the corotation feature becomes narrower and deeper, until 
it reaches a final width. This finite final width is again indicative 
of non-linear behavior, and from Fig.[T2]it is clear that this non- 
linear structure is responsible for the cumulative torque on the 
planet reaching positive values. The ripples on the edges of the 
corotation region in Fig. [12] are of numerical origin and have to 
do with the strong non-linearity in the density gradient. They can 
be reduced by slowly ramping up the mass of the planet (which 
was done for the runs in Fig. [T2b but in the absence of heat dif- 
fusion they never go away completely. 

In Fig. [T3] we show the imaginary part of the to — 2 
Fourier component of the density for the three planet masses 
(black lines). These allow us to make an easy estimate of the 
width of the horseshoe region Xs- Interestingly, if we define 
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Fig. 11. Total torque in units of q''T.{rp)Cl^rp as a function of 
time for 3 different planet masses and an adiabatic equation of 
state with — —1.2 and a = —1/2. The torque is divided by 
to filter out purely linear effects. 



the horseshoe region to extend to the distance where the torque 
is half the minimum, we find that Xg — r^yjiq/ih to within 
a few percent. Also, it can be deduced from Fig. [TTj that the 
time scales for which saturation sets are consistent with libra- 
tion period 87r/(3a;s) set by this value, which one may obtain 
from considering the gravitational and centrifugal equipotentials 
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Fig. 12. Real (dashed lines) and imaginary (solid lines) parts of 
the m = 2 Fourier component of the density perturbation due 
to a g = 3.17 • 10"^ planet after 20 (black) and 10 (red) orbits 
in an adiabatic disc with ^ — —1.2 and a — —1/2. The inset 
shows a close-up on the corotation region. A comparison with 
results obtained from linear theory indicates the form of a linear 
response apart from in the corotation region. 
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Fig. 14. Imaginary parts of the m = 2 Fourier component of the 
density perturbation due to a q = 1.25 • 10^^ planet for different 
resolutions with ^ = —1.2 and a — —1/2, together with a run 
with a steeper radial density profile (a = 1/2, with the same 
value for 7). 
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Fig. 13. Imaginary parts of the m = 2 Fourier component of the 
density perturbation due to a q = 1.25 • 10^^ planet for different 
radial entropy profiles with a — —1/2. For the negative entropy 
slope, we show three different planet masses. 



alone. In this context we comment th a t from a consideration 
of h orseshoe streamlines , iMasset et al.l (|2006|) estimated Xs ^ 
r-p y/ 8rp I $Gps I / ) , where $Gps is the gravitational po- 
tential of the planet evaluated at the point on the horseshoe sepa- 
ratrix closest to the planet (as we use this expression only to indi- 
cate a scaling, we neglect an additional enthalpy contribution to 



the potential). According to this we ex pect Xs — / □ > 



, W2^7rp/(3L) 

for an appropriate length scale L. As iMasset et al.l (12006) point 
out, in practical cases, this should be determined by both the 
softening length and the scale height, H. But note that the soft- 
ening length is generally chosen so that these are closely related. 
For example these authors provide the estimate 



- 1.16rpJqrp/H 



(31) 



This corresponds to L — QA95H which therefore agrees with 
the estimate above which used L = Q.QH — bto within ten per- 
cent. They indicate that their estimate fits locally isothermal sim- 
ulation results at smaller softening lengths and that this should 

be used to give Xs when b ^ 0. 

Incorporating the above ideas, foUowing lMasset et al.l(l2006l) 
from now on we define the horseshoe width through 



'2qrp|$Gps 
3GA/o 



(32) 



where A is expected to be of order unity and can be determined 
from simulations. 

The total torque, divided by q^, depends on the mass of the 
planet (see Fig.fTTli. The reason for this is at least partly through 
the form of the density profile around corotation. For the torque 
to be proportional to , we would expect that the mass associ- 
ated with the feature, or any of its Fourier components, should 
scale as q. From Fig. [13] both the width and the depth of the fea- 
ture are seen to depend on mass but in such a way that the area 
increases more weakly than q^. But note that the small scale 
structure apparent in the density perturbation around corotation 
may well depend strongly on numerics, and therefore also the 
total torque. We will show below that this problem can be over- 
come by adding a small amount of thermal diffusion. 
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Fig. 15. Total torque in units of q^'S{rp)ri^rp as a function of 

time for po cx r~^/^ and q — 1.26 • 10^^. The solid line shows 
an isothermal run, the dotted line an adiabatic run. 
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Fig. 16. Total torque in units of q-'Y,{rp)il,f^rp as a function of 
time on a (7 = 3.17 • 10^^ planet for different values of the heat 
diffusion coefficient D^Ana disc with a = —1/2 and ^ = —1.2. 



In Fig. [14] we show the corotation feature for three different 
resolutions. It appears that the corotation feature gets stronger 
with increasing resolution, and it is not clear if the torque has 
converged for our highest resolution. A common feature for all 
runs is the positive total torque. The total torque is also remark- 
ably robust to variations in the flux limiter, showing differences 
of less that 5% betwee n the min mod limiter and the soft super- 
bee limiter (see Paardekooper & Mellem a 2006b, for their defi- 
nition). From Fig.[T4]it is clear that there are strong features in 
the density that have a width of one grid cell for all resolutions. 
At the two highest resolutions, we resolve all relevant scales: the 
pressure scale height, the Hill radius of the planet and the width 
of the horseshoe region. The appearance of these small scale fea- 
tures even at these high resolutions suggests that the torque is 
limited by numerical diffusion. We will return to this issue be- 
lo w, where we di s cuss m odels with explicit heat diffusion. 

iMasset et al.l (l2006l) noted non-linear behavior in the 
barotropic case above a certain mass. Our analysis suggests that 
in the non-barotropic case, for sufficiently low thermal diffusiv- 
ity, and for sufficiently long integration times, non-linear effects 
are always dominant (see Eqs. ( |28l l and (|29]l). In addition, the 
discussion of the linear problem indicates that, as the state vari- 
ables (but not their gradients) converge as e 0, non-linear 
effects should be much milder in the barotropic case. This is 
confirmed by Fig. |9] where the non-linear evolution takes the 
torque on the planet to positive values only in the non barotropic 
case. 

In Fig. [T3] we show the results for different entropy gradi- 
ents, with the planet mass ratio being fixed at q — 1.26 • 10^'"'. 
The run with no entropy gradient resembles the Lorentz profile 
of the linear case (see the right panel of Fig. |2]l, confirming that 
it is the entropy gradient that drives the non-linearity. A positive 
entropy gradient changes the sign of the corotation feature, mak- 
ing the corotation torque on the planet negative. The run with no 
entropy gradient is intermediate between the positive and neg- 
ative entropy gradient runs and indicates a linear dependence 



of this torque on ^. However, because we measure total torques 
arising from a combination of physical processes, all of which 
vary when the disc structure is modified by altering the entropy 
gradient, this is not conclusive, because entropy-related torque 
contributions may not always be reliably identified. 

The run for which we show the corotation feature in Fig. 
[14] behaves differently when a is changed from —1/2 consid- 
ered so far to a = 1/2. In Fig. [15] we show the evolution of 
the total torque, comparing an isothermal run with an adiabatic 
run with ^ = —0.8. For this disc, the total torque never be- 
comes positive,a indicating a strong dependence of the migra- 
tion rate on the density profile of the unperturbed disc. In fact, 
there are three effects conspiring to make the torque negative in 
the a = 1/2 case. First of all, the wave or Lindblad torques 
are expected to be almost the same being only slightly weaker 
(see Tanaka et al. 2002) . Second, the positive contribution from 
the linear barotropic corotation torque, that would apply to the 
locally isothermal case if the temperature variation across the 
corotation region may be neglected, is expected to be weaker 
by a fact or of two because the radial vortensity gradient is less 
strong ( Goldreich & Tremainelll9'79l) . Not only will the linear 
torque be more negative, there is also no expected strong non- 
linear boost for the barotropic corotation torque in these circum- 
stances. This is clear from the solid line in Fig. [15] which does 
not show the strong rise that was apparent in Figs.|9]and[TT] For 
this example the final torque is also in reasonable agreement with 
that obtained from the linear analysis which is also indicative of 
less important non-linear effects at corotation. Finally the radial 
entropy gradient is less negative than in the case with a = —1/2, 
which reduces the positive contribution of the non-barotropic 
corotation torque. The result is a negative torque on the planet, 
and therefore inward migration, albeit at a slower rate due to the 
positive contribution of the entropy gradient to the torque. 

' Note that in a cylindrical disc the vortensity gradient is zero for 
Po oc r"^/^ (see Eq. fTTll). 
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Fig. 17. Real (dashed lines) and imaginary (solid lines) parts of 
the m = 2 Fourier component of the density perturbation due to 
a q = 3.17 • 10^^ planet after 40 orbits, for different values of 
the heat diffusion coefficient D^, in a disc with a — —1/2 and 
^ = —1.2. The inset shows a close-up on the corotation region. 
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Fig. 18. Total torque in units of q^T^{rp)n'^rp as a function of 
time for po cx r~^/^ and two different planet masses. For the 
q — 3.17 • 10^^ planet, we have used — 7 ■ 10^^, for the 
g = 1.26 • 10"^ planet we have used = 5.6 • 10"^. 



6.3. Effect of heat diffusion 

In this Section we add a small thermal diffusivity into the model, 
characterized by the value of D^. We expect that for a high 
enough value of D^, in the vicinity of the corotation radius, 
all non-linear behavior due a radial entropy gradient disappears 



and the model should behave as if it was locally isothermal. In 
Fig. [16] we show the time dependence of the total torque on a 
q — 3.17 • 10^^ planet for different values of with ^ — —1.2 
and a = —1/2. It is clear that heat diffusion decreases the 
strength of the non-linear part of the torque. However, for times 
less than 5 orbits, the behavior is independent of Dc- This is 
consistent with the solution being in the linear regime during 
this phase (see Sect.[3]above). 

For Dc = 7 ■ 10^^ the torque during the later non-linear 
phase is still positive, but less so than in the adiabatic case 
{Dc = 0). Thus this value of Dc is enough to affect the non- 
linear solution, which is consistent with the estimates made in 
Sect. |3] (see Eq. (|28]l). Increasing Dc by a factor of 10 gives an 
even weaker positive torque on the planet. The fact that after 
5 orbits there is a constant shift in torque between this case and 
the isothermal case indicates that all non-linear effects due to the 
entropy gradient have disappeared. The difference between these 
cases is simply due to different Lindblad or wave torques, which 
is supported by the fact that after 5 orbits the torque ratio is about 
equal to 7. Together with the mildly non-linear barotropic part 
of the corotation torque (which is not affected by heat diffusion, 
but which may still involve significant density gradients in the 
coorbital region, see Sect. [8]l this is enough to switch the sign of 
the torque to positive values. 

To illustrate the behavior of the entropy-related torque, we 
show the rn ~ 2 Fourier component of the density perturbation 
in Fig.[T7]for the cases Dc = 0, Dc = 7-10"^ and £>„ = 7-10-^ 
It is clear that the effects of heat diffusion are localized around 
corotation, where the temperature gradients are strongest. The 
close-up around corotation shows a smooth Lorentz -profile for 
the most diffusive run, resembling the linear case in a qualitative 
way (see the right panel of Fig. |2]i. Note that Eq. ( l28T l predicts 
that non-linear effects should be present for Dc < 3.4 • 10^^, 
while from Fig.[T6]we see that for a diffusion coefficient that is 
twice as high non-linear effects associated with the entropy gra- 
dient are still present. Bearing in mind that Eq. (|28] | gives only a 
rough estimate, and that non-linear effects may indeed set in for 
a somewhat lower values of Dc, the agreement is reasonable. 

Even in the less diffusive case the sharp features at corotation 
that can be seen in Fig.[T2]for the Dc — case are smoothed out. 
In fact, for Dc = 7 ■ 10^^ we find the same density response for 
a simulation with half the resolution, and also the total torques 
agree to within 5 %. Therefore, the torque is numerically well- 
defined, unlike the case with Dc — 0. If we now compare dif- 
ferent masses with a minimum amount of diffusion, enough to 
smooth out the sharp features but at the same time keep the non- 
linear behavior, we find that the total torque scales almost as q^ . 
Note that in order to have the same amount of non-linearity in 
the problem, we need to use different values of Dc for differ- 
ent masses (see Eq. (l28Tl). The results are shown in Fig. [TS] for 
a = 1/2 and ^ — —0.8, where we have rescaled the time for 
the low-mass planet by a factor of 2. This removes the g^^/^ 
scaling for the time to set up the torque. Interestingly, the torque 
plateau in these cases is reduced by a factor of 2 compared to 
the case with Dc = 0, being negative and about 1/2 of the linear 
value. This analysis supports the idea that the departure from the 
scaling of the torque with q^ is due to non-linear effects in the 
corotation region. 

7. A simple non-linear model 

Because linear theory always breaks down for sufficiently small 
Dc, it is necessary to develop a non-linear picture to understand 
the results of the numerical simulations presented above in Sect. 
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|§1 In this section, we discuss a very simple non-linear model that 
can give at least qualitative insight into the torque behavior. It 
is a non-linear model because a finite width of the corotation 
region, in which the gas travels on horseshoe orbits, is adopted. 
Such a region is not present in linear theory. The model only 
accounts for the breakdown of linear theory and the build-up of 
a non-linear contribution to the corotation torque in the early 
stages after a planet is inserted into the disc. The later expected 
saturation of the torque is not considered. 

A description of the way the model works is as follows. 
Viewed in the frame rotating with the planet, material exterior 
to the planet approaches the azimuth of the planet along approx- 
imately linear trajectories. It executes a horseshoe turn, subse- 
quently moving away from the planet when it will eventually 
will have received a radial displacement equal to twice the ini- 
tial distance from the orbital radius of the planet (see the stream- 
lines plotted in Fig.[T9]l. If the planet were absent material would 
shear past without a turn. Thus the effect of the horseshoe turns 
is to produce physical changes in the material that correspond 
to a radial displacement equal to the width of the horseshoe re- 
gion after it has encountered the planet. The horseshoe region is 
separated from the remainder of the disc by a separatrix. 

When the material conserves entropy as it moves, such a dis- 
placement produces a density change. For example, if the disc 
entropy per unit mass decreases outwards, when material is dis- 
placed inwards, the density will increase. Here we shall assume 
the displacement occurs such that fluid near separatrices in the 
horseshoe region maintains pressure balance with its surround- 
ings and by extension that the Eulerian change in the pressure 
may be neglected. In this context we comment that it is found 
both from linear theory and non-linear simulations that the rela- 
tive density changes always greatly exceed the relative pressure 
changes (see above and below). We note that the dynamics we 
have described are clearly seen in the simulations. In Fig. [19] we 
show a close up of the corotation region near the planet for a 4 
M® planet. The density structure indicates an increase inside the 
horseshoe region that is bounded by the inner separatrix. This is 
also indicated by the form of the streamlines superposed on the 
plot. Note that these are slightly distorted by the grid resolu- 
tion. We comment that the change in density that occurs after a 
horseshoe turn is expected to result in a slight asymmetry in the 
positions of the horseshoe legs on either side of the corotation 
radius away from the planet. But for the relative density changes 
considered here (typically of order one percent) this is a very 
minor effect. 

The horseshoe turns commence immediately after a planet 
is inserted into a disc and produce a torque that builds up over 
time in addition to that expected from linear perturbation the- 
ory. Departures from linear theory are thus produced early on. 
As regions close to the planet are the most important in estab- 
lishing the torque, a plateau is reached fairly quickly when these 
regions attain a quasi-steady state. This lasts for the order of a 
libration period after which material recirculates to the location 
of the planet and the torque is expected to saturate. The model 
described here accounts for the attainment of the torque plateau 
but not the final saturation. 

To describe the model quantitatively, we use coordinates 
{x,9) that are centered on the planet: x = (r — rp)/rp, 9 = 
(fi — (fip. As a fluid element turns to move from one side of 
the planet to the other, once it has moved there is ultimately a 
jump from x to —x; after this the streamlines are along lines of 
constant x. Here we neglect minor effects that may arise from 
the slight offset of the centre of the horseshoe region from the 
planet's orbital radius due to the background pressure gradient. 
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Fig. 19. Density, multiplied by r^/^, after 10 orbits of a g = 
1.26-10^''^ planet embedded in a disc with po cx r^^/^, T oz r^^ 
and ^ ~ —0.85. Approximate locations of the streamlines are 
indicated by the black lines. 



For simplicity, we concentrate on the part of the corotation 
torque that is due to the entropy gradient only. In general one 
also expects an additional contribution if the gradient of specific 
vorticity is non-zero. Because most of the torque arises close to 
the planet we expect the major part of it to be established before 
a full libration period and that this could be determined locally. 
Below we estimate the maximum torque arising from the initial 
horseshoe turns as described above. 

Because entropy is conserved along streamlines and we as- 
sume the change Eulerian in pressure is zero, we can write for 
the perturbed density (we focus on a jump from the horseshoe 
leg exterior to the planet that has 6* > to the corresponding leg 
interior to the planet, located in the upper segment of Fig. \l9[ : 

p' = -2xipo, (33) 
7 

where — d\og{P / p^) / d\ogr . The torque on the disc can be 
found by integrating over the disc domain that has participated 
in a horseshoe turn: 

r = -2Hrl j p'^^dOdx. (34) 

The integral over 9 is made simple by adopting as integration 
limits 9 = Q and 9 large enough such that the force due to the 
planet is small, which leaves us with: 

r = 4i/por-p-GMp / ° J_^ xdx, (35) 
7 Jo + ¥ 
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where we have approximated the planet potential by: 



Gp 



Thus we arrive at: 



1 rpb + ^(rp6)2 + 1-2 



(36) 



(37) 



The torque on the planet has the opposite sign and may be writ- 
ten in the form 



^^^s 2y 3q2 



79 (''p^ + \Jxs +rlb'^ 



(38) 



where we have introduced a factor of 2 to account the additional 
contribution arising from the horseshoe legs on the other side 
of the planet. Thus we assume a symmetry such that these be- 
have in a corresponding manner but with a density deficit rather 
than increase produced near the separatrix after a turn, so also 
producing a positive torque. 

We comment here that the above torque depends on the soft- 
ening parameter, b through evaluation of the planet potential 
at the point on the separatrix closest to the planet, here taken 
to be the actual planet location. We may accordingly replace 
it with the magnitude of this potential denoted by |$Gps| = 
GMp/ (rpb). Then the torque may be written in the form 



-fqGMpTp (l 



£il£Gp=lf 



1 



(39) 



In this form the expression for the torque may be used in the 
situation where the point on the separatrix closest to the planet is 
separated or offset from the planet position by an angular extent 
^s- From Eq. ( |36] | we see that this is equivalent to replacing b 

to substitute for Xs- The 



by ^/W+9^. Now we use Eq. 
expression for the torque then becomes 



Tp 



84 
3 7 



1 + A A 



)2£ /rp 
3 V < 



GM„ 



(40) 



where A is the dimensionless parameter defining the horseshoe 
width that is expected to be of order unity that was introduced 
above (see Eq. (|32]|). 

Equation ( |38] ) predicts a torque that is proportional to the 
entropy gradient of the unperturbed disc, and is positive for neg- 
ative entropy gradients. The torque that is reached is comparable 
in magnitude to the isotherm al wave, or Lindbla d, torque for the 
expected value of & w H/r (iTanaka et al.ll2002l) . The torque, as 
given by Eq. (|38] l, does not diverge for 6 — > 0. With any separa- 
tion or offset of the type described above being neglected, setting 
6 = in Eq. ( [38] l gives 



^P = 



-^gSprpl^p. 



(41) 



Note that in this limit the only stream parameter that enters this 
expression is the horseshoe width Xs- The torque is then in- 
creased by a factor Xs/ {rpb 



+ r^b^) compared to the case 

b 0. Using Eq. dSTT l to estimate Xg, this amounts to a factor of 
4 when b = Q.QH/r for q = 1.26 x 10"^ 



However , its us e in this limit is inappropriate. Following 
iMasset et al.l (l2006l) . we note that when vertical stratification is 
considered it would be expected to produce an effective soft- 
ening parameter b = zjvp for material moving at height z. 
Adopting this and performing a vertical average of Eq. (l38l l for 
the above example would indicate b ^ Q.'aH/r, being compara- 
ble to values adopted here. 



7. 1. Time development of the non-linear torque 

We studied the time dependence of the disc response as a result 
of the insertion of a planet of 4 into a disc with po oc r^'^/'^ 
and T oc . The model was adiabatic with De = and 
^ = —0.85, in t his case corresponding to 7 = 1 .1. Note that, 
as mentioned in iPaardekooper & Mellemal (l2008l) . there is no 
smooth transition to the locally isothermal torque for 7^1, 
because in the latter case entropy is not conserved along stream- 
lines (if there is a radial temperature gradient). In fact, for a neg- 
ative density slope, a lower value of 7 tends to make the entropy 
gradient more negative, and therefore the associated corotation 
torque more positive. 

Unlike for the other simulations considered in this paper, in 
order to study the time development of the non-linear torque in 
a way corresponding to the discussion of the simple model, the 
planet was initially immersed in the disc with its full mass. Thus 
there was no period of slow build up for these cases. From the 
streamlines illustrated in Fig. [T9] and discussed above in Sect.]?] 
the inner and outer separatrices for a model with b — 0.03, away 
from the planet, lie on circles centered on the central object with 
dimensionless radii 0.982 and 1.012 respectively. 

In order to study the development of the high-density region 
leading the planet which is responsible for the non-linear torque, 
we plot the density as a function of azimuthal angle, (p, as viewed 
along the radius r ~ 0.99 after 5, 10 and 20 orbits in the right 
panel of Fig. |20] The planet is located at (p = n and for higher 
values of ip, this radial location lies in the required high density 
region. The left panel of the same figure shows the total torque as 
a function of time (black line). A small thermal diffusivity, = 
10^^ was used to regularize the solution around corotation. For 
comparison, we also show the total torque for a simulation with 
Dc = (green line). The torque in this case rises to a slightly 
higher value at t — 10, followed by a slight decay. This effect 
is due to the sudden introduction of the planet into the disc, and 
it can be reduced by introducing a small thermal diffusivity, as 
shown in Fig.|20l or by slowly ramping up the planet's mass. 

We see that as time progresses an approximately uniform 
density increase by a factor of about 1 + 5 x 10^^ is established. 
This is entirely consistent with the factor 6p/p = 1.55x/r, with 
X being the distance to the coorbital radius, expected from the 
discussion above. The uniform density increase is also consis- 
tent with a uniform pressure equal to the unperturbed value be- 
ing attained at the fixed radius. 

Let us now check that the time development is consistent 
with fluid elements moving from one side of the horseshoe re- 
gion to the other Assuming this to be the case, taking account 
only of Keplerian shear so neglecting any time required to cross 
close to the planet, it would take 5 orbits to set up the high den- 
sity region out to ip — 3.4, corresponding to 5 disc scale heights 
from the planet, 10 orbits to set it up out to (p = 3.64 and 20 
orbits out to ip = AAA. However, Fig. |20l shows that the density 
profile is not fully established after 5 orbits, but that after 10 and 
20 orbits the profile is indeed established out to the expected lo- 
cations. This indicates that the evolution is relatively slow only 
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Fig. 20. Left panel: total torque, in units of q'^Y,{rp)Q'^rp, on a g = 1.26- 10 ^ planet embedded in a disc with po (x r T cx r ^ 
and ^ = —0.85. A small thermal diffusion coefficient Dc = 10^^ was used to regularize the solution near corotation. Results for 
two values of the softening parameter b are shown in black and red, together with a simulation with b — 0.03 and Dc ^ (green 
line). Right panel: azimuthal variation of the density at r = 0.99 at three different times for b ~ 0.03 and = 10^^. 



in the regions close to the planet. This could be related to passage 
close to the separatrices where low velocities are expected. 

Note too that half of the total contribution to the corotation 
torque we have considered here is attained after only 5 orbits 
which is much less than the libration period of about 75 orbits 
which is again consistent with most of the contribution coming 
from close to the planet. The magnitude of the non-linear con- 
tribution to the torque predicted by Eq. dJSl l is approximately 
1000 in units of q'^T,{rp)nl_r^, which agrees perfectly with the 
difference between the minimum ( or also approximately the lin- 
ear torque) and the maximum torque in the left panel of Fig. |20l 
This suggests that there is no additional torque cut-off close to 
the planet, since one was not included in the non-linear model 
presented here. From Fig.[T9]can be seen that indeed the corota- 
tion feature extends to about a smoothing length from the planet 
in azimuth. Note too that although the associated relative den- 
sity change is only a few percent, on account of the small radial 
scale the relative change to the density gradient is at least of or- 
der unity. Additional non-linearity not represented in the simple 
model may be present and possibly be responsible for the varia- 
tions of the ultimate ratio of torque to mass ratio squared seen in 
Fig-El The fact that these ratios are affected by small diffusion 
coefficients in the manner illustrated in Fig. [TS] indicates local- 
ized effects in the coorbital zone are responsible. A possibility 
is that because density jumps across the horseshoe separatrix are 
higher for larger masses, because the horseshoe region is wider, 
non-linear effects are also relatively stronger. 

Also shown in the left panel of Fig.|20]is the total torque for 
a run with a smaller value of the smoothing length b = 0.02. 
From Eq. ( |40l i, with A = 1, we expect the non-linear torque to 
be a factor 2 stronger in this case, and this is entirely consistent 
with Fig. |20] We point out that from this trend the total torque 
would be expected to be always negative for b > H/ r. 
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Fig. 21. Long-term evolution of the total torque, in units of 
g^I](rp)ri|,rp, on a q = 1.26 • 10^^ planet in a disc with 

Po oc r^^/^, T (X r^^ and ^ = —0.85, for three different ther- 
mal diffusivities. 



8. Long - t e r m evolu ti o n 

We now turn to the question of the behavior of the torque over 
longer time scales, focusing on a disc with no vortensity gra- 
dient. We initialize the temperature T cx r~^, and again take 
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Fig. 22. Density, multiplied by r'^/^, for three models for q — 1.26 -10 ^, pa oc r ■^Z^, T cx r ^ and £^ — —0.85 after 160 orbits, 
focusing on the horseshoe turn leading the planet. Left panel: — 10^^, i' ^ 0. Middle panel: ~ 10^^, = 0. Right panel: 



^ = —0.85. As we need to follow the evolution of the torque for 
several libration cycles, which, for a g = 1.26 • 10^^, amounts 
to several hundred orbital periods, we lowered the resolution by 
a factor 2 in all directions to keep the problem computationally 
tractable. As mentioned in Sect. 16.31 for a finite value for the 
thermal diffusion coefficient the torques agreed very well for the 
first 30 orbits. We focus on a (7 = 1.26 • 10^^ planet throughout, 
and, to keep matters as simple as possible, do not ramp up the 
mass of the planet over the first 3 dynamical time scales as was 
done in the previous sections. Again, for a finite thermal diffu- 
sion coefficient, this does not influence the torques at all, while 
for Dc = differences remain well within 10 %. 

In order for a diffusive process to prevent the corotation 
torque from saturating, it is necessary that the diffusion time 
scale across the horseshoe region is smaller than the libration 
time scale, or: 



(42) 



Using our estimate Eq. ([32]) for Xs, which has been shown to 
be a good approximation, we find for q = 1.26 • 10^^ and 
b = 0.03 that > 3/8{xs/rpf « 3 • 10"^ to aflow ther- 
mal diffusion the possibility of restoring the entropy gradient 
within one libration cycle. On the other hand, according to Fig. 
[TSl for Dc — 5.6 • 10^6 the disc response is still sufficiently non- 
linear that the torque is changed considerably. Therefore, there 
exists a range of Dc for which there may be a possibility of a 
sustained positive torque contribution from corotation effects re- 
lated to the entropy gradient. However, thermal diffusion only 
affects the temperature gradient directly, and it may be the case 
that a consequent adjustment of the density in the horseshoe re- 
gion allows it to settle towards a state of zero torque contribution. 

In Fig. im we show the long-term evolution of the total 
torque for three values of Dc- The case of zero diffusion shows 
the clearest libration cycles, diminishing in amplitude until the 
system is expected to settle into a state of vanishing entropy gra- 
dient in the mean and therefore zero corot ation torque. T his be- 
havior is very similar to that described by IWardI (l2007l) for the 



barotropic case with a radial vortensity gradient. Interestingly, 
both diffusive runs approach the same value of the torque, which 
indicates that indeed the density is redistributed in such a way 
that the corotation torque vanishes. 

The final negative torque of magnitude ~ 1200 compares 
with the linear value (based on the orig inal density d i stribu tion) 
of ^ 700. For comparison the result of iTanaka et all (12002) ob- 
tained by use of Eq. ( |30] | for the original density distribution 
gives a magnitude of ~ 900. As indicated above this was in 
much better agreement with the torque obtained for the locally 
isothermal case with the same original density distribution. The 
most likely reason for the greater discrepancy in this case is that 
the original density gradient in the coorbital region and its neigh- 
borhood has been modified as indicated in Fig. |22]in an asym- 
metric manner. 

In the final state, the corotation region contains an asymmet- 
ric density structure. This can be seen in the two leftmost panels 
of Fig. |22] where we show the density after 160 orbits for the 
cases with low and high diffusion. In the leftmost panel, it is 
clear that the low density ridge created in the horseshoe turn be- 
low the planet returns on the other side. For the higher value of 
Dc, a similar situation arises, but since diffusion for this case is 
actually strong enough to affect the density jump directly, the ef- 
fect is harder to spot. From Fig. [21] however, we see that in this 
case approximately the same total torque is attained. 

All of this indicates that in order to sustain the positive 
torque, we need to restore the initial density profile as well as 
the initial temperature profile in the coorbital region. This can 
be done by including a small kinematic viscosity. In Fig.|23] we 
show the torque evolution for the same models as in Fig. [21] 
but now including a viscosity v — 10^^ ''p^^p- The mass flow 
induced by this small value of is not strong enough to affect 
the results; in essence the viscosity only acts on the sharp den- 
sity features near corotation. We see that for a small viscosity in 
combination with a small thermal diffusivity the positive torque 
can be sustained. The magnitude strongly depends on the mag- 
nitude of the diffusion coefficients, as is illustrated by the dashed 
line in Fig.|23] A simulation with no heat diffusion shows simi- 



S.-J. Paardekooper & J.C.B. Papaloizou: Disc planet interactions in non-barotropic discs 



19 



O 



-500 



-1000 - 




150 
t (orbits) 



300 



Fig. 23. Long-term evolution of the total torque, in units of 
g^S(rp)r2|^rp, on a g = 1.26 • 10^^ planet in a disc with 

po oc r^'^/^, T oc r^^ and 7 = 1.1, for three different ther- 
mal diffusivities. All models have a small kinematic viscosity of 
u = 10-6 rlQp. 



lar behavior as the runs depicted in Fig.lJT] but the final value of 
the torque agrees better with linear theory, presumably because 
of the smoother density profile. This indicates that indeed we 
need to restore both the original temperature gradient (through 
heat diffusion) and smooth the density gradient (though a small 
viscosity) in the coorbital region in order to sustain the positive 
torque. 

These results indicate that the entropy-related corotation 
torque alone is indeed able to reverse the sign of the total torque 
over longer time scales than the libration period. Although the 
torque in Fig.|23]is only slightly positive at later times, note that 
this model represents a case for which the entropy-torque alone 
has to overcome the full Lindblad torque. A higher value of 7, 
which reduces the Lindblad torque, and a flatter density profile, 
which gives rise to a positive corotation contribution due to the 
vortensity gradient, may help to push the torque to higher posi- 
tive values. On the other hand increasing 7 may act to reduce the 
entropy gradient. The investigation of how the vortensity and en- 
tropy corotation torques operate together will be the subject of a 
future study. 



9. Discussion 

In this paper, we have considered the torques induced by disc- 
planet interaction when low-mass planets are introduced into 
protoplanetary discs with an initial radial entropy gradient and 
either solved the energy equation or adopted an adiabatic or lo- 
cally isothermal equation of state. We first considered the linear 
theory of the interaction. As the response is singular it needs 
to be regularized. This was done by adopting the well-known 
Landau prescription, together with a small thermal diffusiv- 
ity in some cases. This allows migration torques on the planet 
to be calculated, but the density response diverges as dissipa- 
tive effects are reduced to zero implying that non-linear effects 



will ultimately be important. We estimated that non-linear ef- 
fects would be significant for dimensionless thermal diffusivities 
De < 10^ '^-6-' and that in such cases the linear theory would 
only be valid for some finite time after insertion of the planet. 
This was confirmed by the non-linear simulations. It was found 
from these that a non-linear density structure developed around 
corotation at early times, that in some cases with a negative en- 
tropy gradient, led to a reversal of the sign of the torque leading 
to outward migration. This positive torque survives for a horse- 
shoe libration period after which saturation occurs and the torque 
again becomes negative. 

We presented a simple non-linear model to account for the 
non-linear coorbital density structures as being produced after 
the horseshoe turns undergone by material while conserving its 
entropy in Sect. [T] Numerically we found that in cases with no 
diffusivity, the coorbital density structure showed variations on 
the grid scale and therefore in order to get numerical conver- 
gence at limited resolution, a small non-zero thermal diffusivity 
needs to be included. Increasing the thermal diffusivity brings 
the migration torque closer to the linear prediction although 
some nonlinearity associated with the torque connected with the 
specific vorticity gradient may remain. 

The non-linear and temporary nature of the physical effect 
that causes the torque to become positive makes it inappropri- 
ate to attempt to deduce an appli cable linear torque forrnula, a s 
was done in the barotropic case dGoldreich & TremaindligTOh . 
All that can be said is that this non-linear part of the torque de- 
pends on the entropy gradient. We have found that, for a disc 
with H/r — 0.05, a density gradient less steep than po (x 
results in a positive torque on the planet. However, this depends 
also on the adopted value of the smoothing length, the value of 
7 and the amount of thermal diffusion. A smaller value for b will 
enhance such non-linear effects, while thermal diffusion acts in 
the opposite direction. 

Our simple non-linear model of Sect. Q does a good 
job in qualitatively describing the results of the numeri- 
cal simulations. Because t he thre e-dimensional simulations of 
IPaardekooper & Mellemal (l2006ah show the same slow rise of 
the torque, which is indicative of non-linear effects, it is likely 
that a similar model applies in 3D. 

The strong dependence of low-mass planet migration on the 
radial entropy gradient in the gas disc raises two important ques- 
tions. First of all: do realistic, self-consistent disc models with 
detailed radiative transfer allow for negative entropy gradients? 
Second: is the thermal diffusivity low enough to support non- 
linearity at corotation? The second question was already an - 
swered by the simulations of lPaardekooper & Mellemal (l2006ah . 
which showed that, at least in the inner parts of protoplanetary 
discs, radiative transport does not provide enough diffusion to 
restore linearity. As opacity in general decreases with radius, the 
situation will be different in the outer parts of the disc. 

The first question is more difficult to answer Detailed mod- 
els of D'Alessio et al. (1998, 1999, 2001) show a rich variety 
in temperature and density profiles, with indeed several cases 
with strongly negative entropy gradients, especially near opacity 
jumps. These may be preferential r adii for halting planet migra- 
tion, a conclusion also reached by Menou & Goodman (2004) 
for different reasons. Although it is beyond the scope of this pa- 
per to investigate the occurrence of negative entropy gradients 
in detail, it is important to realize that they do appear in detailed 
models. 



In the inviscid limit, the corotation region is a closed system, 
and can therefore not continue to exert a torque on the planet in- 
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definitely. The first signs of saturation of the torque are already 
visible in, for example. Fig. |9] and the complete process is il- 
lustrated in Fig. |2T] In order to sustain this torque for longer 
times, the physical conditions in the coorbital region need to be 
regenerated. Heat diffusion can act to restore the original en- 
tropy gradient, but we have found that the form of the surface 
density gradient also needs to be smoothed locally, which can be 
achieved by including a small kinematic viscosity. 

At this point we stress the sensitivity of torque behavior de- 
rived from effects at corotation to small values of therm al diffu- 
sivity and viscosity. This has already been indicated by i Massetl 
(1200 2) who found that small values of viscosity acting in con- 
junction with effects at corotation could have significant effects 
on the migration of planets but in a larger mass range than that 
considered here. 
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10. Summary and conclusion 

In this paper, we have studied the effect of including the energy 
equation in the theory of disc planet interactions undergone by 
low-mass planets. We have shown that the wave torque on the 
planet is reduced by a factor 7 compared to locally isothermal 
discs. It is the corotation torque, however, that is the most strik- 
ingly different from the isothermal case. 

We have found that migration behavior depends critically on 
the radial entropy gradient in the disc. The linear estimate of 
the corotation torque depends on this, but the entropy gradient 
is also responsible for a departure from linearity, giving rise to 
a strong corotation torque that is positive for a negative entropy 
gradient. For a strong enough entropy gradient these corotation 
torques can dominate over the wave torques. This can be under- 
stood in a qualitative way by considering simplified horseshoe 
dynamics. Thermal diffusion reduces non-linearity, and for high 
enough values of the thermal diffusion coefficient Dc we recover 
modified linear wave torques and a mildly non-linear barotropic 
corotation torque. 

Although for a planet on a fixed orbit and with no viscos- 
ity, the positive torques saturate after a libration period, we have 
shown that, for a small thermal diffusivity and viscosity the 
entropy-related corotation torque offers the possibility of a sus- 
tained reversal of the direction or stalling of migration of low- 
mass protoplanets. Of course extensive future investigations that 
take into account three dimensional effects and allow the planet 
to migrate are required in order to evaluate it. 
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